Mon. Not. R. Astron. Soc. 000, 000-000 (0000) Printed 1 August 2008 (MN KTbX style file v2.2) 



Effects of gravitational-wave recoil on the dynamics 
and growth of supermassive black holes 



Laura Blecha and Abraham Loeb 

Harvard University, Department of Astronomy, 60 Garden St., Cambridge, MA 02138, USA 
Email: lblecha@cfa.harvard.edu, aloeb@cfa.harvard.edu 



1 August 2008 



00 
O 

o 

< 



6 



(N 
> 
O 



ABSTRACT 

Simulations of binary black hole mergers indicate that asymmetrical gravitational wave (GW) 
emission can cause black holes to recoil at speeds up to thousands of km s^^. These GW 
recoil events can dramatically affect the coevolution of recoiling supermassive black holes 
(SMBHs) and their host galaxies. However, theoretical studies of SMBH-galaxy evolution 
almost always assume a stationary central black hole. In light of the numerical results on 
GW recoil velocities, we relax that assumption here and consider the consequences of recoil 
for SMBH evolution. We follow the trajectories of SMBHs ejected in a smooth background 
potential that includes both a stellar bulge and a multi-component gaseous disk. In addition, 
we calculate the accretion rate onto the SMBH as a function of time using a hybrid prescription 
of viscous (a-disk) and Bondi accretion. We find that recoil kicks between 100 km s^^ and 
the escape speed cause SMBHs to wander through the galaxy and halo for ^ 10^ — 10^ 
yr before settling back to the galactic center. However, the mass accreted during this time 
is roughly constant at ~ 10% of the initial mass, independent of the recoil velocity. This 
indicates that while large recoils may disrupt active galactic nuclei feedback processes, recoil 
itself is an effective means of regulating SMBH growth. Recoiling SMBHs may be observable 
as spatially or kinematically offset quasars, but finding such systems could be challenging, 
because the largest offsets correspond to the shortest quasar lifetimes. 

Key words: black hole physics - gravitational waves - accretion, accretion disks - galaxies: 
kinematics and dynamics - galaxies: evolution - galaxies: active 
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1 INTRODUCTION 

Supermassive black holes (SMBHs, > 10® Mq) have long been 
considered the probable central engines for quasars ( |Lynden-Bell| 
1 1969 ). Observations now indicate that they are ubiquitous in the 
centers of both active and quiescent galaxie s (e.g. |Kormendy| 
|& Richstone] ([T995|; [Richstone eF^ ([T998|; [Ferrarese & Ford| 
( |2005) ). Considerable evidence links SMBHs closely to dynami- 
cal properties of their host galaxies, indicating that SMBHs and 
galaxies may have evolved concurrently. Measured SMBH masses 
are closely correlated with the luminosity, and hence the mass, of 
their host galactic bulges (Kormendy & Rich stone||1995[ Magor-| 
Irian et al.| 1998) [Marconi & Hunt|2003i|Merritt'&'Fe rrarese'2001|l. 
Also, a strong correlation has been observed between SMBH mass 
and stellar velocity dispersion in the b ulge: Mbh oc at jGebhardt| 
let al.|2000[ [Ferrarese & Merritt|2000"| [Tremaine et al.|2()02a^ . Re^ 
cent surveys have found that a similar correlation may extend down 
to dwarf galaxies if their compact nuclei are interpreted as the ana- 
log of SMBHs in massive galaxies l |Ferrarese et al.|2006[|Wehner| 
|& Harris|20()6l|Barth et al.l2005l l. 

These findings indicate that SMBHs offer a wealth of infor- 
mation about the formation and evolution of structure in our uni- 



verse. In this context, much theoretical work has been devoted to 
studying the details of SMBH accretion and how it relates to galaxy 
evolution (e.g.,|Silk & Rees|1998|[Wyithe & Loeb|2003||Di Mat^ 
|teo et al.||2005| [Hopkins et al.||2006| l. However, virtually all such 
studies have one common assumption: that the black hole (BH) re- 
mains stationary at the center of the galaxy. This is a reasonable 
assumption in light of the observed ubiquity of central compact ob- 
jects in galaxies, although the observed sample contains only local 
galaxies, so in principle few constraints exist on the dynamics of 
SMBHs over the lifetimes of galaxies. However, if two SMBHs co- 
alesce after a major galaxy merger, gravitational-wave (GW) recoil 
can impart a large kick to the resultant merged BH ( |Peres|[T962[ 
|Bekenstein|1973[ [Fitchett & Detweiler|1984^ . When the merging 
BHs have unequal masses or spins, asymmetrical emission of grav- 
itational waves creates a net linear momentum flux at merger, caus- 
ing the merged BH to recoil in the opposite direction. This effect 
has long been known as an interesting relativistic phenomenon, but 
its importance as an astrophysical phenomenon remained very un- 
certain until recently; the best estimates from approximate meth- 
ods indicated that kick velocities were most likely to be small 



< few X 100 km s i Blanchet et al. 



[mar 2006| l. Results from simuiations using tull numerical relativity 
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have demonstrated that if high BH spins are included, kicks up to 
4000 km are possible with certain mass and spin configura- 
tions jCampanelli et al.]|2007a|b) . These velocities easily exceed 
the escape speed of any known galaxy, so the ramifications of large 
recoil kicks would clearly be dramatic. Not only will the position 
of the BH in the galaxy vary, but also the rate of fueling and the 
amount and type of active galactic nucleus (AGN) feedback gener- 
ated. For example, according to the Bondi accretion model, which 
is commonly used in simulations, the accretion rate is inversely 
proportional to the velocity of the accretor cubed ( [Bondi & Hoyle| 
[T9^ . 

The velocity of recoil kicks depends strongly on the spin pa- 
rameters, spin alignment, and mass ratio of the merging BHs. The 
distributions of these quantities in the population of merging galax- 
ies are highly uncertain. |Bogdanovic et aLlpOO V) have suggested 
that torques in a circumbinary gas disk may align the BH spins 
with the disk orbital plane, causing the recoil distribution to be 
skewed toward in-plane kicks (and consequently low kick veloc- 
iti es).[Schn ittman & Buonanno] j2007[ l, [Campanelli et"aL]j2007a| l, 
& |Baker et al.| (i2008 ) have each calculated probability distributions 
of recoil kick velocities assuming randomly distributed BH mass 
ratios and spin orientations. We allow for both aligned and random 
spins by considering a variety of kick inclinations. Observing sig- 
natures of GW recoil in the future would be one way to constrain 
these distributions. 

A common result of early SMBH merger simulations was 
the so-called "final parsec problem", whereby inspiralling SMBHs 
"stall" at a characteristic radius of about a parsec after hollow- 
ing out the host galaxy's core and before gravitational wave (GW) 
emission becomes dominant (e.g., [Milosavljevic & Merritt 200 1| 
|Yu|2002^ . This has been referred to as a "problem" because until 
very recently (see [Rodriguez et al.|[2006) , no direct evidence for 
any such unmerged binary SMBHs existed, let alone a large pop- 
ulation of them. Numerous viable solutions to this problem have 
been introduced, for example by adding a slightly triaxial potential 
or a nonnegligible gas fraction (e.g., Gerhard & Binneyll985 [[Yu 
2002 ; Berczik et al. 2006] [Escala et al.||20041 [Gou ld & Miralda- 



Escude|1997jl. It is certainly possible, however, that in "dry" merg- 



ers of galaxies containing small gas fractions, the SMBH inspiral 
time would exceed the time for another major merger to occur. This 
could result in formation of a triple SMBH system, which may im- 
part a significant kick to one of the BHs (e.g., [Hut & Rees[1992[[Xu[ 
[& Ostriker[1994[ l. Recently, the first fully relativistic simulations of 
close triple-BH encounters demonstrated that substantial kicks can 
indeed o ccur jCampanelli et aL|2008] l. Using non-relativistic sim- 
ulations, |Ho^ma^^^oeb[([2007]rconducted a statistical analysis 
of the triple-SMBH problem; they found that most three-body en- 
counters cause the binary BH to merge without ejecting the third 
BH, and that ejection of the third BH (which may also be followed 
by the merger of the other two) occurs in less than half of all en- 
counters. It should be noted that if triple-SMBH systems do form 
in galaxies and if the most likely outcome is an SMBH merger, 
then GW recoil kicks could be at least as likely to result from this 
configuration as dynamical kicks to the third BH. 

[Madau & Quataert| p004| l and [Loeb[ j2007) have considered 
some of the observational signatures we might expect from these 
sources, such as spatially and kinematically offset quasars. These 
could be either SMBHs that are ejected from their host galaxies and 
carry an accretion disk along, or wandering SMBHs on bound tra- 
jectories that may significantly disrupt galactic structure. Because 
no known quasars have been conclusively identified as such, we 
can assume the fraction of offset quasars is relatively small. Simi- 



larly, [^lonten]([2007| considered the limits placed on ejected BHs 
by the Mbh — o"» relation. She concluded that ejection of SMBHs 
could not have been a very common occurrence throughout most 
of the history of the universe, or we would observe considerably 
more scatter in BH-galaxy relations. However, Volonteri points out 
that in large halos at high redshift, mergers were more frequent and 
recoil kicks, which are independent of absolute mass, could more 
easily have ejected BHs from the relatively small host galaxies. 

[Gualandris & Merritt| ( [2008^ studied GW recoil with a set of 
N-body simulations following the motion of a kicked SMBH in 
a stellar potential. They find that BHs on bound orbits have long 
oscillation timescales, up to ~ 1 Gyr, including a phase of low- 
amplitude oscillations in the stellar core. A large stellar core is 
scoured out by these oscillations, in some cases larger than the 
cores that binary SMBHs are expected to produce. 'Guala ndris &[ 
Merritt ( 2008 1 also discuss the possible observational signatures of 
such oscillating BHs, such as displaced AGN and offset jets. 

In a recent paper, [Komreich & LoveIaceH2008y have focused 
on a different aspect of GW recoil, namely its effect on galaxy mor- 
phology. Using GADGET-2, they also find that BHs on bound or- 
bits have long oscillation timescales and conclude that recoil events 
can produce significant asymmetry in galactic disks. We note, how- 
ever, that by definition, any galaxy containing a recoiling SMBH 
will have recently undergone a merger, which will also produce 
asymmetric morphology. 

Magain et al.^(2005( ) discovered a quasar with no apparent host 
galaxy located near an Ultra-Luminous Infrared Galaxy (ULIRG), 
which had likely undergone a recent merger A number of papers 
discussed whether this might be an ejected SMBH jHaehnelt et al.[ 
p006; Merritt et al. 2006, Hoffman & Loeb 2006'). This does not 
appear to be the case for this system, however. Among other fac- 
tors, the quasar has a significant narrow-line region indicating the 
presence of more gas than could be ejected with the SMBH, and re- 
cently, direct evidence for a host galaxy has been observed jPeam] 
|etaL|20 07). 

While this paper was undergoing revision, [Komossa et al.[ 
l [2008| announced that the quasar SDSSJ0927-I-2943 is a candidate 
recoiling SMBH with a velocity of 2650 km s~^. They have ob- 
served an offset between the broad emission lines and some of the 
narrow emission lines in this unusual quasar's spectrum, which they 
suggest can be explained if the broad line region is gravitationally 
bound to a recoiling BH and the narrow line region is left behind 
in the galactic disk. 2650 km s^^ is close to the maximum possi- 
ble recoil kick velocity, and the actual velocity may be even higher 
if there is a velocity component perpendicular to the line of sight. 
Thus, this proposed recoil event is certainly not a typical event. 

[Bonning et al.[ ( [2007) conducted a targeted search for GW re- 
coil events. They examined a quasar population for Doppler shifts 
in spectral lines that would indicate a moving source and found 
a null result. However, this method applies only to high-velocity 
BHs that are active quasar sources. Additionally, because a BH on 
a bound trajectory spends a large fraction of its time at turnaround, 
such objects are statistically more likely to be observed with low 
velocities. 

In this paper, we focus on the dynamical aspect of SMBH ac- 
cretion in the context of large kicks from GW recoil. We emphasize 
that the details of BH accretion occur on scales below the current 
resolution of state-of-the-art simulations of galaxy mergers, and 
that the physics of turbulent viscosity driving this accretion cannot 
be simulated easily. Thus, a semi-analytic approach is an appropri- 
ate first step toward understanding accretion onto recoiling BHs. 
We use a set of semi-analytical models to follow the trajectories 
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of SMBHs ejected with varying kick velocities in tlie potential of 
a galaxy containing both gas and stellar components; we also es- 
timate the accretion rate for the moving SMBHs. Our model and 
choice of parameters are outlined in § |2] The results of these cal- 
culations for our fiducial model are given in § [3ll while §|3]2]& 
1 3. 3 1 detail the results of varying the SMBH mass and galaxy gas 
fraction. In § |3.4[ we consider the implications of our results for 
observational signatures of GW recoil. § |3.5| contains comments on 
the sensitivity of our results to our choice of parameters. Finally, in 
§[4] we summarize and discuss these results in the broader context 
of SMBH-galaxy coevolution. 
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2 MODEL 



We construct a galaxy model consisting of two distinct compo- 
nents: a spherical stellar bulge and a gaseous disk. Note that in real- 
ity, the stellar and gas distributions may be highly asymmetric fol- 
lowing a major merger. Choosing symmetric distributions is there- 
fore a simplification necessary for our semi-analytical approach, 
although if the SMBH merger timescale is long, the galaxy may 
in fact be significantly relaxed when the recoil occurs. The dark 
matter (DM) halo is ignored for the purpose of following SMBH 
orbits, because the halo profile will be fairly flat on the scale of 
the baryonic galaxy. Also, on large scales the halo is expected to 
be significantly triaxial jjing & Suto|2002[ [Bullock et al.|2001) , so 
even though the DM halo would decrease the SMBH orbital decay 
time, it would also increase the chance of scattering the BH onto a 
less-centrophilic orbit, such that the BH would wander in the halo 
instead of returning to the galaxy center. The details of this scatter- 
ing are highly dependent on the shape of the triaxial potential, so 
we ignore this effect and consider it a source of uncertainty in our 
high-velocity recoil runs. 

2.1 Stellar bulge 

The stellar bulge is assumed to follow an 77-model density profile 
( [Tremaine et al.|2002b with 77 = 2.5. We choose this profile that 
is slightly flatter than a Hemquist] \ 1 990| > profile (77 = 2) under the 
assumption that a galaxy merger remnant will have a flattened core 
due to stellar ejections during SMBH binary inspiral. The total stel- 
lar bulge mass is estimated from the SMBH mass using the relation 



of |Marconi & Hunt|p003) , 



M* = 4.06 X 10'"M 



/ Mbh 



(1) 



The scale radius, a, for the r;-model profile is related to the SMBH 
mass in terms of the half-mass radius, ri/2, 

a = ri/2(2^/'' - 1). (2) 

An approximate relation for the ratio of effective (half-light) radius 
to half-mass radius is given by |Dehnen| ( |1993[ Eq. (17)) 

Rc 



ri/2 



0.7549 - 0.00439(3 - rf) + 0.00322(3 - ■qY 



-0.00182(3 - r?)'^ ± 0.0007. (3) 

The effective radius is i?o = GA-h/3a^, and the one-dimensional 
velocity dispersion, a, , is calculated from the Mbh — cr, relation 
( [Tremaine et al.|2002a^ , 

0.249 



(J, = 220 km 



Mb 



2 X IO^Mq 



Figure 1. Schematic face-on diagram of our model for the gaseous galactic 
disk (not to scale). Zone I is the standard fShakura & Sunyaev]jl973) «-disk, 
marked with horizontal hash-marks. The solid circle denotes the radius rqi . 
Zone II (shown in dark gray) is the inner (BH-dominated) portion of the 
self-gravitating accretion disk, while Zone III (shown in light gray) is the 
outer (disk-dominated) portion of the accretion disk. The transition between 
the two, rtrans, is marked with the long-dashed circle. The BH radius of 
influence, ring, is marked with the dot-dashed circle. Beyond this radius is 
the exponential disk (Zone IV), which has no defined outer boundary. 



Thus, the scaling of the stellar model radii with BH mass is (i?e, 
ri/2, a) oc Mg^j^. 

Stellar dynamical friction on the ejected SMBH is modeled 
using the [Chandrasekhar] ( | 1 943| l formula. Assuming a Maxwellian 
velocity distribution, the drag force on a BH moving with velocity 
vbh through a background density is: 



fdf = -I{M) X — '-^-^ , 

<jt «BH 



where 



TT 



(5) 



(6) 



(4) 



and the Mach number A4 = vbh/o"*. We use 3.1 for the Coulomb 
logarithm ln(A) ( jEscala et al.|2004[[GuSandris & Merritt|2008[ >. 



2.2 Gaseous disk component 

In addition to a stellar bulge, most galaxies will contain a significant 
gaseous component. This is especially true of merger remnants, as 
mergers are more frequent between galaxies at high redshift, where 
gas fractions are generally higher. We assume that by the time a 
SMBH binary has merged in the galaxy center, the gas has cooled 
into a disk. For our fiducial model, we assume the disk mass is 
determined by the gas fraction /gas = Mdisk/M* = 0.5, such that 
the galaxy is relatively gas-rich. 



2.2.1 Disk structure 

Analytical models of disk structure generally fall into two distinct 
categories: SMBH-dominated, viscosity-driven accretion disks 
(i.e., Qf-disks) with Mdisk ^ Mbh, and self-gravitating, star- 
forming galactic disks with Mdisk S> Mbh. The intermediate 
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regime, Mdisk ~ Mbh, is less straightforward to model analyti- 
cally and is seldom resolved by hydrodynamical simulations. This 
regime also defines the SMBH radius of influence, rinfl, through 



the relation Mdisk (< nnfl) = 2Mbh. Goodman (2003 & Good 



[man & Tan| ( |2004} discuss the issue of self-gravitating accretion 
disks, though they focus on scales of < pc where the disk is still 
nearly Keplerian (as in our Zone II, defined below). Thompson' 
|et al.|p005| conduct a more detailed analysis of this problem all 
the way to the BH radius of influence using a semi-analytical model 
of the disk structure that includes dust; this is a more involved ap- 
proach than is necessary for the scope of our paper 

In addition to seeking a simpler disk model, we address the 
more general question of whether a purely analytic disk solution 
can be constructed for the entire region between the self-gravity ra- 
dius and the BH radius of influence. Solutions have been derived by 
other authors for various limits, but to our knowledge this is the first 
attempt to extrapolate these limits to model the entire disk. Our ap- 
proach is to construct an analytic model of the disk structure at all 
radii, from the innermost stable circular orbit around the BH to the 
outer reaches of the galaxy. Because this range incorporates vastly 
different physical environments that cannot be modeled by a sin- 
gle set of simple equations, we define several regimes separated by 
critical radii. Fig.[T|shows the four concentric "zones" we have de- 
fined in this "onion-skin" model. Extending the disk continuously 
from galactic scales down to the innermost radii is important for 
calculating the accretion rate and ejected mass for a wide range of 
recoil velocities. However, we note that BHs kicked out of the disk 
plane spend a small fraction of their wandering time in the inner 
disk regions, such that the dynamics of these BHs will generally be 
dominated by the outermost disk region (Zone IV) and the stellar 
distribution. Thus, our dynamical results are largely independent of 
the details of the inner disk model described here. 

The innermost zone. Zone I, is a standard a-disk l |Shakura"&l 
|Sunyaev|T9 73 ). The outer edge of this region, rqi, is defined as the 
radius at which the disk becomes unstable according to the |Toomrel 
( |I964) parameter. 



Q 



ttGE' 



(7) 



where Cs is the gas sound speed, 



„2 



+ r{dn'^/dr) 



the epicyclic frequency defined in terms of the rotation angular 
frequency, fi, and E is the disk surface density. Q < 1 corre- 
sponds an unstable (i.e., self-gravitating) disk. Following |Good-| 
[man & Tan| ( |2004} , we assume radiation pressure dominates gas 
pressure at r < rqi, such that the sound speed is simply given by 
Cg = AgbT* /icp, where ctb is the Stefan-Boltzmann constant. 
Although radiation pressure dominates in this region, we assume 
the viscosity scales only with the gas pressure; otherwise the disk 
may be unstable. From these assumptions one can derive the radius 
rqi ( |Goodman|2003j [Goodman & Tan|2004| : 



rqi 1.5 X 10 pc 



a/0.3 



1/3 



/e 



Erad/O.l 



1/6 



VKcs/ 



1/2 



IO8M0 



1/2 



dependence of rqi on BH mass is weak; thus, it is reasonable to 
calculate rqi assuming that Zone I of the disk is radiation-pressure 
dominated in general (see Table[TJ. Likewise, for r > rqi the disk 
is assumed to be gas-pressure dominated. 

The self-gravitating region beyond rqi will cool and undergo 
star formation, which will heat the disk via feedback processes. Gas 
compression and shocks may also form due to gravitational insta- 
bilities created when Q < 1. It has frequently been argued that 
these feedback and cooling processes will balance, leaving the disk 
in a marginally stable state with Q ~ 1 (e.g., |Hohl|I97I|[Paczynski] 
|I978 [ Bertin 1 997). This is su pported by results from various sim- 
ulations (e.g., |Gammie|200r||Lodato & Rice|2004^ . Accordingly, 
we constrain the disk beyond rqi to have Q = 1, which is effec- 
tively an equation balancing heating and cooling processes in the 
gas. We derive our set of equations using this constraint in place of 
the equation for balance between radiation and viscous dissipation 
of orbital energy (see, e.g., Bertin 1997; Bertin & Lodato I999j 
|Goodmanl[2003| ). We retain the assumption of a steady accretion 
rate via viscous dissipation of angular momentum. 



(8) 



Here /i is the mean molecular weight of the gas, /Edd ~ 
Ma/ M-Edd is the ratio of the a-disk accretion rate to the Edding- 
ton rate, trad is the radiative efficiency, and the opacity k is scaled 
to the electron-scattering value, = 0.4. We use k = Kcs in 
our calculations. For our fiducial parameters, the radius at which 
the radiation and gas pressures are equal is rpoq ~ rqi, and the 



Ma = SvTi^E = STTQfCs^E = const, 



(9) 



where the viscosity v = ac^h is defined according to the a-disk 
prescription, h is the disk scale height, and the dimensionless con- 
stant a is a free parameter, usually assumed to be in the range 
~ 0.01 — 0.1. |Goodman| p003] l suggests that for self-gravitating, 
Q = 1 disks, this value may be as high as 0.3. This is the fiducial 
value we use in our disk calculations. Consumption of gas via star 
formation in the disk is neglected, such that M remains constant, 
although the energy input from star formation is included indirectly 
by setting Q = 1. We also use the thin-disk approximation, as well 
the relation pg = E/2^ for the gas density. 

We construct the intermediate disk regime as follows. Just 
beyond rqi, the BH still overwhelmingly dominates the gravita- 
tional potential. Likewise, at the BH radius of influence, rinfl, the 
disk becomes dominated by its own gravitational potential. We use 
these two limits to define two zones (Zones II & III) in the region 
tqi < 1' < ''infl: in Zone II the disk contribution to the poten- 
tial is neglected, and in Zone III the BH term is neglected. The 
transition radius, rtrans, is defined where the two density profiles 
intersect (En = Em), thus creating a smooth transition between 
the two regimes. In the thin-disk, BH-dominated regime (Zone II), 
Cs = /ifisH, where J^bh ~ GMbh/t^ ■ The relation for Cs, along 
with Eqs. |7| & l|9]l, can be solved to give expressions for E, Cs and 
h: 



f^BH 
TV 



Ma 



1/3 



Cs.II 



/ GQMa 

\ 3a 



1/3 



hll = il-n 



f GQMa 

\ 3a 



1/3 



(10) 



(11) 



(12) 



For Zone III, Cs is approximated by using only the disk poten- 
tial term to solve for vertical hydrostatic equilibrium and the radial 
Poisson equation, giving Cg — 27rGE/i. 0{r) similarly is taken to 
depend only on the disk mass, Mdisk: n{rf = GMdisk(< r)/r\ 
This is the form of f2(r)^ for a spherical distribution; using this 
expression for a disk is initially a simplifying assumption. How- 
ever, we find a self-consistent solution for a disk with a 1/r surface 
density profile ~ i.e., a Mestel disk (Mest el|1963) . In this special 
case, the disk's rotation curve exactly equals that of an equivalent 
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spherical distribution, so is not an approximation. The general 
equation for the viscous mass accretion rate in a non-Keplerian disk 
is: 



d Inr 



= const. 



(13) 



Eqs. ^ & along with the sound speed equation, yield the disk 
quantities for Zone III: 



Ca.III = 



4 



M 



1/3 



1/3 



hi 



8 



(14) 



(15) 



(16) 



This results in an angular vel ocity i ^ jr)^ ~ 
{csjii/hui)^ (hiu/r). In a recent review, 'Lodato| 1 2008 also 
derives the set of equations for what we call Zone III and defines 
the radius of transition to a Keplerian disk. Note that his numerical 
factors differ slightly from ours due to a different definition of the 
disk scale height. 

We model the region beyond rinfl (Zone IV) with an exponen- 
tial profile, Eiv ~ Eoe"'"'^ ' where the scale radius rdisk is set 
by matching Siii(nnfl) = Eiv(rinfl)- So is obtained by integrat- 
ing E(r) from rinfl and equating this to the gas mass beyond rinfl: 



Eo = 



A/disk (> nnfl) 



'-disk J 



ntl/'-dii 



(17) 



To calculate the disk gravitational potential at points outside 
the disk plane, we again use the approximation of spherical sym- 
metry for the gravitational potential, namely, d^/dR = GM{< 
R)/B?. However, a perfectly spherically symmetric potential is 
clearly unphysical given the asphericity of the gas disk. To break 
this symmetry, we introduce a small elliptical perturbation to the 
potential by writing the radial coordinate as _R = r"^ + {z/qY, 
where q < 1 is the axial ratio of an equipotential ellipsoid. This pa- 
rameter is not intended to accurately reproduce the ellipticity pro- 
file of the galactic potential at all radii - in fact, in this approxima- 
tion the ellipticity is constant at all radii. Instead, it merely allows 
some precession of the orbits in the r — z plane, as one expects 
due to torques introduced by the aspherical mass distribution. As a 
fiducial value in our runs, we choose q = 0.99 (i.e, a 1% distortion 
of the axial ratio). We find that the largest effects due to this pertur- 
bation occur for near-escape orbits with ~ Ucsc, where the BH 
reaches the largest radii. This further validates the choice of small 
q, since in reality the galactic potential will tend toward sphericity 
at large radii. See § |3.5| for further discussion on the sensitivity of 
our results to this parameter. 

When the SMBH is ejected from the center of the galaxy, it 
carries with it a portion of the gaseous disk out to a radius rcj . This 
radius can be estimated simply as r-cj = GA/bh / (I'k + ) ; namely, 
the gas orbiting the BH with VoA > «k remains bound to the BH 
as it leaves the galaxy center jLoeb|2007) . We thus assume the total 
mass ejected is Afcj = Afan + Afdisk.cj. Another consideration is 
that numerous studies have found that in hydrodynamical simula- 
tions a binary BH with moderate mass ratio (such as would pro- 
duce large recoils) creates torques on the surrounding gas disk that 
carve out a circumbinary gap about twice the size of the semimajor 
axis (Liu et al. 2003 , Escala et al. 2005 ; Milosavljevic & Phinney 
|2005, MacFadyen & Milosavljevic, 2006, ,Hayasaki et al.n2007| l. 



As the binary approaches merger, the gravitational-wave merger 
timescale, tgw, decreases rapidly, and the BH orbit decays faster 
than the hole's viscous timescale to refill, ivisc. The final gap size 
at merger is determined by the point at which tgw << tvisc, which 
for our models is ~ 100 rs ( |Milosavljevic & Phinney|2005| rs is 
the Schwarzschild radius). The disk mass within this radius is neg- 
ligible, so we can ignore the circumbinary gap in our simulations. 

We also assume that the ejected disk leaves a hole in the galac- 
tic gas disk of radius in our model. The hole will eventually 
refill; we assume refilling will take place on the viscous timescale 
tvisc ~ r'^ /v = /acsh, though we note that Milosavljevic & 



|Phinney| ( |2005[ l have suggested that the hole should refill faster 
than the viscous timescale. We calculate this timescale in each of 
our runs, based on the environment at the hole radius rej . For BHs 
kicked into the disk plane with low velocity, this timescale will 
be > 10* yr, which in this case is generally longer than tfln, the 
time for the BH to settle back to the center. For out-of-plane kicks, 
^visc > ifln for low ujj, whcu Tej is large. The velocity V]^ where 
iviac = ifln depends on the galaxy model and kick inclination, and 
is ~ 300 — 400 km s~^ for our fiducial model. We neglect the re- 
filling of the central hole in these cases. However, for larger kick 
velocities, the hole is very small (rej oc and makes little dif- 

ference to the overall dynamics and accretion of the moving BH. 



lU- 

I 



2.2.2 Gas dynamical friction 

When the BH + ejected disk are moving through the galactic disk, 
gas dynamical friction acts on them in addition to stellar dynam- 
ical friction. 'OstrikerJ ( |1999| l derived an analytical formula for a 
body moving through a gaseous medium that has the same form 
as Chandrasehkar's formula, but with a different function I(A^) 
and with M = ubh/cs. This formula provides a fairly accurate 
estimate o f gas drag force s and has been used by many other au 
thors (e.g., 'Narayan 2000; 'Karas & Subr"200l' 'Kim et al."2005' 
However, Escala et al. ( 2004 1, as well as Sanchez-Salcedo & Bran 
|denburg| ( |2001^ , have found that Ostriker's formula overestimates 
the dynamical friction somewhat for slightly supersonic velocities, 
~ 1. To get around this problem, we adopt the approach of |Es-| 
|cala et aL| ||2004 ) and choose a simple parametric form for I(A^) 
that approximates Ostriker's formula for a given value of ln(A) 
but lowers the drag force in the region M ~ 1 — 2. Assuming 
ln(A) = 3.1 as a reasonable nominal value (cf. |Lin & Tremaine] 
|1983[ [Cora et al.|1997^ , we use Chandrasekhar's formula (Eq. [5]l 
for the gas dynamical friction, but with ln(A) = 4.7 for M > 0.8 
andln(A) = 1.5forX < 0.8. 

We must also consider that the kinetic energy removed from 
the SMBH by dynamical friction is added to the gas disk; this 
energy input is not negligible. Detailed treatment of this process 
would require a full hydrodynamical approach, but ignoring it com- 
pletely yields unphysical results. We therefore modify our prescrip- 
tion for gas dynamical friction to account for the heating of the disk 
as follows. We assume the disk is "puffed up" to a modified scale 
height h' ~ r, such that the gas density is reduced by a factor h/r 
and the sound speed Cs is increased by a factor (h/r)^^, where h 
is the initial scale height. We furthermore assume, due to the large 
initial energy input as the BH moves through the cold disk, that this 
heating occurs effectively instantaneously, and we always use the 
modified quantities in the dynamical friction calculation. 
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2.2.3 SMBH mass accretion rate 

In addition to tracking the trajectory of the SMBH, we also make 
an estimate of the accretion rate at each timestep. Nominally, the 
accretion rate is set to Ma, a free parameter for r < rinfl (disk 
zones I - III). We define Ma — /Edd x AfEdd, where 

AttGMbh 



MEdd = 



(18) 



f^cs C trad 

and we adopt a value of e^ad = 0.1 for the radiative efficiency. We 
set /Edd ~ 1 for the gas-rich (/gas ~ 0.5) models, such that the 
accretion is Eddington-limited. For consistency we keep the same 
ratio /Edd = 2/gas for the /gas = 0.1 model (such that /Edd = 0.2 
in this case). If the SMBH is out of the plane of the disk, its ac- 
cretion rate is Ma for < t < taccr, where facer is the accretion 
timescale Afdisk.cj /Ma ■ After this time, we assume the ejected disk 
has been consumed, and the accretion rate is set to zero. When the 
SMBH is in the gas disk, either as it leaves the center or on sub- 
sequent passages through the disk ("disk crossings"), the accretion 
rate may change. Again, att < taccr, we can nominally assume that 
the BH carries with it a disk of its own that feeds the BH at a rate 
Ma. However, the BH may pass through regions of the galactic gas 
disk with a much larger scale height than in its own disk, meaning 
that the BH may accrete mass from all directions, not just its disk 
plane ~ i.e., mass in the galactic disk may be swept up by the BH. 
The accretion from above and below is more accurately modeled 
by Bondi-Hoyle accretion ( |Bondi & Hoyle|1944| l: 

47r(GAf)Vg,oo 



Mb = 



(19) 



where Vroi(r) = [ t;BH,r(r)f, (hbh, i;circ) 0, wbh,z(j')z ], 
^t-ij;. = iij-ij;. (r + Tej ) is the circular velocity at the outer edge of the 
BH disk, and Pg,oo, Cs,oo are also evaluated at this radius. This can 
be rewritten in terms of the Bondi radius re = 2Gi\'fcj/(i'roi + Cs) 
as 



Mb = vr/9g,oo \/ Ujci + ci.cx) rlfg, 



(20) 



where we have added the geometric factor /goom = miri(l, /i/re) 
to account for the reduced cross-sectional area of accretion if the 
Bondi radius is larger than the physical height of the disk. Because 
the BH also carries with it a disk from which it is accreting at a rate 
Ma, the accretion cannot be modeled by either the Bondi-Hoyle 
or a prescriptions alone; we therefore introduce a hybrid prescrip- 
tion. We first assume that as long as t < facer, the BH accretion 
rate will not fall below Ala while passing through the galactic disk, 
because it accretes at this rate from its bound disk. Bondi-Hoyle 
accretion alone decreases with BH velocity as l/tjjljj, but it also 
varies with environment, which depends on radius. It is unclear a 
priori whether the accretion rate will increase or decrease in this 
case. We therefore calculate the accretion rate as follows. MB(r,\) 
for a moving object is related to A/b(0), the Bondi rate for a sta- 
tionary, central object by 



.(0) 



MB(r,v) = f X Mb{0), 

3 



\MbhJ 



Pg,oo(0) 

/geom 
/geom(O) 



(21) 



For a cold disk rather than an idealized infinite, homogeneous back- 
ground, the low gaseous sound speed in the disk gives a Bondi ra- 
dius much larger than the BH radius of influence when the BH is 



stationary, while in reality the two should be about equal. To avoid 
calculating unphysical accretion rates much greater than AfEdd, 
we define a fiducial sound speed Cs,oo = {2GM/rintiY^^ , where 
M = Mbh for Cs,oo(0) and A/ = Afej for Cs.oo- These are the 
sound speeds used in the equation above. Because we know that 
the BH will nominally accrete at the a-model rate instead of the 
Bondi rate, we adopt the factor / defined in Eq. l[2TJ, which de- 
pends on Cs,oo, i^roi, Pg,tx), Afcj, and /gcom, and normalize it to 
Ma. We then have the accretion rate M — /amp x Ma, where 
/amp = max(/, 1) prevents M from falling below Ma. As men- 
tioned above, this applies only when t < facer; otherwise, the BH 
is naked when it passes through the galactic disk and its accretion 



rate is just AfB(r, v) as given in Eq. 1 21 



From this approach, and also adding the requirement that the 
accretion rate be Eddington-limited, we can define five distinct ac- 
cretion regimes for the BH: 

• |2:bh| > h,t < facer : M = Ma 

• I^BhI > /!.,f > facer : A/ = 

• IjZbhI < h,t < facer, /amp = 1 : M = Ma 

• |2Bh1 <h,t< facer, /amp > 1 : A/ =min(/amp A/q, , AfEdd) 

• |zbh| < h,t > facer : M =min( A/b (r , v) , Affidd ) 

In practice, we find that the amplification factor /amp is al- 
most never greater than unity; slight amplification occurs in our 
dry-merger model for the (low-velocity) in-plane kicks. In all other 
cases the v^'^ dependence in the Bondi-Hoyle formula dominates 
the dependence on the changing gaseous environment, such that 
any contributions to the mass accretion from above and below the 
BH-disk plane are negligible. In other words, a negligible amount 
of mass is swept up by the BH. We can therefore choose /Edd = 1 
for our fiducial model, such that Ala ~ A/sdd, and disregard the 
amplification factor. When f > facer, however, M is determined 
solely by the Bondi-Hoyle formula and can have more interesting 
behavior as seen in Figs.[3]|4] &|9] 

2.3 Calculation of SMBH trajectories 

We integrate the trajectories of SMBHs under the influence of grav- 
ity and dynamical friction for a range of kick velocities and angles. 
The integration is stopped when one of four criteria is met: 

• The BH escapes from the galaxy. The actual escape velocity is 
slightly higher than that calculated from the gravitational potential 
alone, due to the drag from dynamical friction. We consider the BH 
"escaped" when it exceeds the nominal escape velocity, \/— 2<1?, 
and is still moving away from the galaxy at a distance of 5 rdisk- 

• The BH settles back in the center of the galaxy, with velocity 
vbh < O.Olcr,, where ct* = a, (Re) is the stellar velocity disper- 
sion calculated from the A/bh — o"» relation. 

• Sufficient energy is injected into the gas via dynamical fric- 
tion to effectively unbind the disk. The binding energy of the gas 
disk is estimated to be ~ 0.5 G Afdisk/(5 rdisk). The cumulative 
energy injected by dynamical friction is the sum of the energy at 
each timestep df, ~ A/cj fdf vbh df. (This limit is not reached in 
any of our simulations.) 

• The integration time reaches a Hubble time. (This limit is not 
reached in any of our simulations.) 



2.4 Fiducial model parameters 

In our fiducial model (Model A), we use a BH mass of 10* M© and 
a gas fraction of /gas = 0.5. Other galactic parameters are scaled 
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Model 




A 


B 


c 


D 


Mbh 


[Mq] 


108 


106 


109 


109 


fgas 




0.5 


0.5 


0.5 


0.1 


M* 


[IQIOMq] 


4.1 


0.034 


45 


45 




[10-2 pc] 


4.0 


0.026 


50 


15 


rqi 


[10-2 pc] 


1.5 


0.15 


4.7 


3.6 


'^trans 


[pc] 


13 


2.8 


28 


82 


''infl 


[pc] 


60 


6.9 


160 


250 


fdisk 


[102 pc] 


7.8 


1.1 


20 


18 


a 


[102 pc] 


7.2 


0.60 


25 


25 


MEdd 


[M0 yr-1] 


2.2 


0.022 


22 


22 


Ma 


[Mq yr-i] 


2.2 


0.022 


22 


4.4 



Table 1. Parameters used in different galaxy models. The free parameters, 
A/bh 2nd /gas. are shown in bold. The other parameters are derived from 
these as defined in the text. 



to these two input parameters as outlined in § |2.1| & [2^ and Table[T| 
summarizes these quantities. We test kick inclinations of i = 0°, 
45°, and 90° relative to the disk plan^ Two merging BHs that have 
no spin or aligned spins oriented along the binary orbital rotation 
axis produce a recoil kick in the binary orbital plane. This very 
specific configuration would be important only if the BH spins are 
preferentially aligned by some mechanism, such as torques from 
the gas disk (Bogda novic et al.|2007^ . In this case we can assume 
that the orbital plane and resultant recoil kick would both be aligned 
with the disk plane. Simulations using full numerical relativity have 
shown t hat the maxim um recoil kick velocity in this case is ~ 200 
km I Gonzalez et al. 2007b I, so we include only data with ?;k < 
200 km s- for in-plane (i = 0°) runs. 

In the absence of empirical evidence for preferentially aligned 
spins, we must also consider that the distribution of merging BH 
spins may be closer to random, producing a random distribution of 



kick inclinations (cf. Schnittman & Buonanno 2007 , Campanelli, 
let al.|2007a| [Baker et al.|2008| l. Thus, we also consider kicks with 
i = 45° & i = 90°. For these out-of-plane kicks, recoil veloci- 
ties co uld be very large ( [Campanelli et al.|2007a[ [Gonzalez et al.| 
|2007 a!| - up to 4000 km s-^ (Campanelli et al.|2007 bV However, 
trajectories are uninteresting for > Vesc, so we calculate trajec- 
tories for kick velocities ranging from 100 km s-^ up to the escape 
velocity, in increments of 20 km . 

2.5 Additional models 

We also test three models in addition to our fiducial model (see 
TablefTJ. We calculate the same set of trajectories with BH masses 
of lO'^nd lO'' M0 (Models B & C, respectively), representing the 

^ All kicks velocities are oriented along the positive x-axis; thus, inclined 
kicks are initially in the x-z plane and in-plane kicks are initially purely on 
the X-axis. 



lower and higher ends of the observed SMBH mass function. 10^ 
M0 BHs are of particular interest because the mergers of these BHs 
would be observable with the LISA gravitational wave detecto^ In 
addition, because models A, B, & C are "gas-rich" models, we also 
test a model with A/bh = 10^ Mq and a gas fraction of /gas ~ 
0.1 (Model D), corresponding to a massive, low-z galaxy that has 
tmdergone a nearly "dry" merger. 



3 RESULTS 

3.1 Fiducial model 

All the results in this subsection pertain to our fiducial model 
(Model A). As explained in § |2j only low velocities ( < 200 km 
s-^) are important for in-plane (i — 0°) kicks. Very little hap- 
pens in these runs; the kick is effectively a small perturbation that 
is quickly damped out in < 10^'* yr. We do not expect these low- 
velocity, in-plane kicks to be easily observable. 

Inclined kicks with low kick velocity are also mostly unin- 
teresting; an example is shown in Fig. |2] Here, the trajectory and 
mass accretion rate are plotted for the i = 45°, Uk = 100 km s^^ 
run. The BH (and its accompanying disk) wander about 20 pc away 
from the galactic center in the x and z directions, settling back to 
the center in < 4 Myr after just one orbital period. The timescale 
for the BH to accrete the ejected disk, taccr, exceeds the wandering 
time; hence, the accretion rate AI = Ma is constant throughout the 
simulation. 

We have plotted an example of a BH with a moderate kick 
velocity (Wk = 440km s"^) in Fig.js] Compared to the low-velocity 
example in Fig. |2] the wandering time of the BH is much longer, 
~ 15 Myr, and and the BH travels much further from the galactic 
center. Note also that the size of the ejected disk, and hence taccr, 
are smaller for higher kick velocities. In this case, taccr ~ 9 Myr 
is marked with a red dot in Fig. [Sj), and corresponds to the point at 
which AI first drops to zero. The subsequent spikes correspond to 
passages through the disk, or disk crossings, where the BH is again 
fed by the gas via Bondi-Hoyle accretion. 

The effective escape speed for i = 45° kicks is between 740- 
760 km s-^. At Hk = 740 kms-^ (Fig. [4|, the behavior is a qualita- 
tively similar but greatly amplified version of the moderate-velocity 
example in Fig. [3] The orbital precession induced by the ellipticity 
of the galactic potential is quite apparent in this case. The BH con- 
sumes its small ejected disk almost immediately and reaches a max- 
imum distance of ~ 14 kpc from the center. During the ~ 1.2 Gyr 
it takes to settle back to the center, it crosses the disk 28 times, pro- 
ducing multiple short bursts of accretion (Fig. |4]5). This is an ide- 
alized scenario, however, because on scales > kpc, the DM halo 
potential will start to dominate that of the baryonic galaxy com- 
ponent. Because the DM distribution is likely to be significantly 
triaxial on these scales, the BH may be scattered onto an orbit that 
takes significantly longer to return to the center of the galaxy, or it 
may wander indefinitely in the halo. (See |Vicari et al.|2007[ for a 
more detailed analysis of this regime.) This means that the number 
of disk crossings and total mass accretion are upper limits, and the 
BH wandering time is a lower limit. The amount of accretion will 
in reality be reduced even further; 1.2 Gyr is well in excess of typ- 
ical starburst timescales (~ 100 Myr, cf. e.g., [Springel et al.|2005[ 
[Marcillac et al.|2006^ , after which less gas will be left to accrete. 
Additionally, this wandering time is roughly comparable to typical 
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Figure 2. (a) SMBH trajectory for Model A run with i = 45°, f k = 100 
km s~^, seen from three orthogonal orientations. Turnaround points are 
numbered in sequence to show the BH trajectory. The j/-axis (top & side 
panels) has been expanded by a factor of 10 to show greater detail, (b) top 
panel: radial distance of SMBH from center vs. time. Maxima are numbered 
to correspond with turnaround points in (a). Bottom panel: Calculated mass 
accretion rate vs. time, normalized to the a-model accretion rate. Note that 
in this low-Ilk example, the accretion rate is constant (M = Ma) for the 
entire run. 



Figure 3. Same data as Fig.|2] but with = 440 km s~^. In Fig.|3^, the 
y-axis (top & side panels) has been expanded by a factor of 100 to show 
greater detail. Turnaround points are numbered in sequence to show the 
BH trajectory; only the first 3 are numbered. The red dot in the top panel 
of Fig. 3p indicates the time taccr = ^^disk,oj /^^q when the ejected disk 
has been consumed by the SMBH. Correspondingly, the accretion rate (bot- 
tom panel) drops to zero here, and is nonzero again only when the SMBH 
crosses the galactic disk. 



timescales between major galaxy mergers, so this may be another 
complicating factor. We note, however, that the upper limit of ~ 1 
Gyr for the wandering timescale in massive galaxies is a similar 
result to that of |Gualandris & Merritt| ( (2008| l, despite the substan- 
tially different physics and parameters that were included in our 
respective models. 

When the BH is kicked directly out of the plane (i = 90°), 
the effective escape velocity is slightly lower than in the i = 45° 



case (~ 660 km s"^ for Model A). The main difference in this case, 
however, is that very little accretion occurs after facer • This is be- 
cause the BH precesses very little about the galaxy center when it 
is ejected directly along the vertical axis; all of its disk crossings 
occur very close to the center where the disk height is small, so the 
disk crossing time is typically too short for any measurable accre- 
tion to occur. 

To understand how GW recoil may affect the mass accre- 
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Figure 4. Same data as Fig.|2] but with tijj = 740 km (i.e., just below 
the escape speed for this inclination from this galaxy model). In Fig.|4^, the 
j/-axis (top & side panels) has been expanded by a factor of 100 to show 
greater detail. Turnaround points are numbered in sequence to show the BH 
trajectory; only the first 9 are numbered. The red dot indicates the timescale 
taccr, as described in Fig.|3] 



tion history of SMBHs, we calculate the total mass accreted by 
the SMBH during each simulation (AM). We are also interested 
in how this compares with the mass that would have been ac- 
creted by a stationary BH of the same mass over the same time 
(AMo = ten/Ma). Fig.jsjshows the SMBH growth as a function 
of kick velocity. AMq, and correspondingly tfi^, increase mono- 
tonically with kick velocity until Vesc is reached. In the middle 
panel where AM is shown, one can see the remarkable result 
that the mass accreted by wandering SMBHs is fairly constant at 
~ 10% /or all kick velocities ^Wcsc- This appears to be a coin- 
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Figure 5. SMBH mass increment for Model A, as a function of kick veloc- 
ity. In all panels, black filled squares with dotted line denote i = 0° kicks 
(low Vk only), blue filled circles with thick line denote i = 45° kicks, and 
red filled triangles with thin line denote kicks with i = 90°. Open cir- 
cles around a point indicate that the BH escaped from the galaxy's baryonic 
component. Top panel: The estimated BH growth, normalized to the initial 
BH mass, if the BH were to remain in the galactic center for the simulation 
time (AA/q = Ma x tfln). Because AMq scales linearly with tgn, the 
latter is shown on the right vertical axis. Middle panel: The estimated BH 
growth, normalized to the inital BH mass, calculated based on the moving 
BH's trajectory. Bottom panel: The ratio of these two quantities. 



cidental balance between the competing factors of the wandering 
time, which increases with Uk, and the time to accrete the ejected 
disk, which decreases with higher «k. § |3.4| gives further details 
about this result. 



3.2 SMBH mass models 

In Model B with a 10*' M© BH, kick velocities as low as ^ 200 km 
s~^ result in escape of the BH from the galaxy. Since the magnitude 
of GW recoil depends on the mass ratio of the progenitor BHs but 
not on their absolute masses, the probability of large recoil kicks 
is the same for smaller SMBHs as it is for the most massive ones. 
Note that AM /Mbh ~ 10% for v < Vcsc in this model as well. For 
V > Vcsc, however, one can see that AM/A/bh gradually declines, 
because taccr is also decreasing (see Fig.|l l[l. 

Fig.jrjshows that the growth for a 1(FM0 BH is qualitatively 
similar to that of a 10* Mq BH, but here the trends are even more 
apparent. AGN feedback or gas depletion via star formation should 
reduce the accretion onto stationary BHs over the long timescales 
shown in the top panel; the projected growth factor of almost 100 
is clearly unphysical. However, the middle panel demonstrates that 
recoil kicks are effective at suppressing the BH growth; the growth 
factor is again remarkably constant at ^ few x 10% regardless of 
kick velocity. 
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Figure 6. SMBH growth for Model B as a function of kick velocity. Nota- 
tion is the same as in Fig.|5] 
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Figure 8. SMBH growth for Model D as a function of kick velocity. Nota- 
tion is the same as in Fig.|5] 
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Figure 7. SMBH growth for Model C as a function of kick velocity. Nota- 
tion is the same as in Fig.|5] 



3.3 Dry merger model 

Model D (a "dry merger") has the same parameters as the lO'' 
model above, but with a gas fraction of 0.1. Fig. [8] shows the BH 
growth for this model, which is again qualitatively similar to Fig. [7] 
As would be expected, the escape velocity is lower for the galaxy 
containing less gas (note that we have defined /gas such that this 
model, with lower /gas, is less massive overall than the others). 
Aside from this, the growth factor for the wandering BH is still 
~ 10%, suggesting again that dynamics, not the details of the gas 




1.5xl0» 



I [yr] 



Figure 9. Same notation as Fig.[2j), here showing data from Model D for 
a trajectory with Dk = 1100 km s~^, i = 45°. In this case, Bondi-Hoyle 
accretion amplifies M above Ma = 0.2ME(jd up to the Eddington limit. 



physics, dominate the determination of the final BH mass when the 
BH is not stationary. 

Fig.|9]shows an example of the BH accretion in Model D, with 
Ala < M^dd- The kick velocity — 1100 km s^^ is just below 
the effective escape velocity for this model, so taccr is very short 
and AI quickly drops from Ma to zero. The accretion rate on all 
subsequent passages through the disk is set by the Bondi rate Mb, 
which varies depending on the velocity and position of the disk 
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crossing. Near the end of the simulation, as the BH slows down, its 
accretion becomes Eddington-limited. 



3.4 Observational signatures 

Because GW recoil has not been unambiguously observed in any 
astrophysical system, an important application of this work is to 
consider the observational signatures that recoiling BHs could pro- 
duce such that their existence could be directly confirmed. An ob- 
servable recoiling BH, i.e. one with significant accretion luminos- 
ity, could be distinguished from a stationary central BH via offsets 
in either physical space or velocity space. A BH with a resolvable 
spatial offset from its host galaxy is statistically more likely to have 
a relatively low velocity, since a large fraction of orbital time is 
spent at turnaround. Likewise, the largest velocity offsets will oc- 
cur either soon after the recoil event or on subsequent passages of 
the BH through the galactic disk, where the spatial offset from the 
galactic center is small. The duty cycle for the former type of event 
depends on the size of the accretion disk that the recoiling BH car- 
ries with it; higher recoil speeds result in smaller ejected disks and 
shorter duty cycles. After this disk is exhausted, the most plausible 
way for the wandering BH to become observable again is by ac- 
creting matter as it makes subsequent passages through the galactic 
disk, possibly generating a quasi-periodic signature such as knot- 
ted or twisted jets. Accordingly, we distinguish between two types 
of "recoiling quasars:" the former, primarily spatially-offset active 
BHs are identified as "off-center quasars," and the latter, primar- 
ily velocity-offset sources are identified as "disk-crossing quasars." 
Note, however, that there is not a one-to-one correspondence be- 
tween quasar velocity offsets and disk-crossings, nor between spa- 
tial offsets and "off-center" quasars. In particular, ejected quasars 
will often only be visible for a short time after ejection, while they 
still carry an accretion disk, so that off-center quasars will often 
have a significant velocity offset as well. Similarly, many "disk- 
crossing" events occur when the BH is settling back to the center 
of the galaxy and slowing down, so their velocity offsets will be 
small. Note also that the term "quasar" is used loosely here to re- 
fer to a source with observable accretion activity; recall that in our 
runs all BHs except those in Model D radiate at their Eddington 
luminosities. 

To quantify the duty cycles for off-center quasars, we plot both 
the simulation time fgn (the time for the BH to settle back to the 
galactic center) and the calculated accretion time for the ejected 
BH disk, taccr, as functions of V]^ (Fig. 1 10^. The bold lines show 
the actual time for which the BH might be visible as an off-center 
quasar. At low kick velocities, the BH settles back to the center 
before consuming its bound disk, and at high kick velocities, the 
BH accretes its disk and then continues to wander as a naked BH. 
The effective time for off-center quasar activity, ~ 10^'^ — 10^ yr, 
varies relatively little with Uk. This balance helps explain why the 
growth of the wandering BHs is also roughly constant with varying 
kick velocity (Figs.|5]|6]|7][8j. However, Fig. |11| shows that taccr 
continues to decrease, albeit slowly, for > 800 km s~^. Thus, 
we expect a lower AM/Mbh for 3> fcsc. 

The simple picture shown in Fig. [TO] provides a plausible ex- 
planation for the accretion behavior seen in our simulations based 
on off-center quasar activity alone, suggesting that accretion dur- 
ing disk crossings is relatively unimportant. Also, for Uk ~ fcsc, 
AM/Mbh ~ 10% regardless of whether the BH escapes or re- 
turns for subsequent disk passages, further indicating that disk 
crossings contribute negligibly to the growth of recoiling BHs. 
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Figure 10. Relevant timescales for off-center quasars. Black line: lifetime 
taccr = A/(jisk,cj/A^Q of the accretion disk ejected with the BH. Blue line 
with filled circles: simulation time tgn for runs with i = 45°. Red line with 
filled triangles: tfj^ for runs with i = 90° . In all cases, the bold portions of 
the line indicate the minimum of (iaccr, tfln) at each i)k ; remaining portions 
are shown with a dotted line. 
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Figure 11. taccr = Mdisk,cj/A/Q is plotted for 100 < Uk < 4000 km 
s^^, for models A, B, C, & D. The solid dot on each curve denotes the 
escape velocity for that model. The bump in curve B occurs when r^j < 
rqi and the disk profile changes. 
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Figure 12. Number of disk crossings versus for eacli model. Blue with 
filled circles indicates kicks with i = 45°, and red with filled triangles 
indicates i = 90°. Upper left: Model A. Lower Left: Model B. Upper 
right: Model C. Lower right: Model D. Only bound trajectories are shown, 
because the escaping BHs do not cross the disk. 

However, a more quantitative analysis demonstrates that this is not 
the complete picture. 

As one would expect, the number of disk crossings (A^'cross) 
generally increases with increasing kick velocity (Fig.|12^. In each 
model, even Model B with the least massive BH, Naoss > 10 
for the highest kick velocities. The kick inclination has little effect 
on the number of disk crossings. However, as mentioned in § |3.1| 
the amount of quasar activity during these crossings is inclination- 
dependent; virtually no accretion occurs when i = 90°, due to the 
small disk height in the inner regions. Additionally, very low kick 
inclinations will produce few disk crossings; we therefore expect 
recoil kicks with moderate inclinations to be the most likely can- 
didates for producing observable disk-crossing quasar activity. We 
accordingly restrict our discussion of disk-crossing quasars to our 
runs with i — 45°. 

Figs. 1 13^ & b compare the contributions of both quasar types 
to the total quasar activity. Fig. |13^ shows the total mass accreted 
for the fiducial model, as well as the mass accreted via off-center 
quasars only and via disk crossings only. We see that as expected, 
the off-center quasar phase dominates the overall growth, although 
at kick velocities 550kms~^ <«k < '^csc, the two contributions 
are comparable. Note that in this same velocity range, the simu- 
lation time is generally greater than the expected timescale for a 
merger-induced starburst, ~ 10* yr. Because accretion during disk 
crossings does not dominate the BH growth, however, a paucity of 
fuel in the post-starburst phase would reduce the total mass accre- 
tion by only a moderate amount. Furthermore, if triggering of the 
starburst happens late in the merger process, the fuel supply may 
last longer. 

Fig.|13|3 is the same as Fig. [To] except that the i = 90° runs 
are excluded, and an extra line is added to denote the total active 
quasar time for each simulation, taccr + icross- Interestingly, the to- 
tal time spent in disk-crossing quasar phases is quite large at high 
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Figure 13. (a) AM/Mbh for Model A, i = 45° . Blue thick line with filled 
circles: total mass accreted during simulation (same as in Fig. |5}. Open 
circles denote simulations in which BH escaped from the galaxy. Green 
thin line with crosses: mass accreted from the ejected accretion disk, as an 
off-center quasar. Magenta dashed line with crosses: mass accreted during 
disk crossings, (b) Same plot as Fig. |10[ but with i = 90° runs omitted and 
a (cyan) line with crosses added denoting the total quasar duty cycle for 
each simulation, i.e, taccr + icross, where tcross is the sum of duty cycles 
for all disk crossings in each run. 



velocities, up to ~ 100 Myr in a few cases. This result has two 
important caveats, however. First, the accretion rate is clearly not 
Eddington-limited for the entirety of these long duty cycles, since 
the mass contribution from this type of accretion is at most compa- 
rable to that of off-center quasars. Evidence for this sub-Eddington 
accretion can be seen in Fig.|4|). Correspondingly, the luminosities 
of disk-crossing quasar events may often be relatively low and thus 
harder to detect than off-center quasar activity. The second caveat 
is that the longest duty cycles for disk crossings are those in which 
the BH has the lowest velocity, i.e., the final passages before the BH 
settles back to the galaxy center. This is a tradeoff in observability: 
the disk-crossing quasars that would be most easily distinguished 
observationally from stationary quasars are those with the highest- 
velocity disk passages and, by definition, the shortest duty cycles. 

Fig.[T4]demonstrates this point more quantitatively; the distri- 
butions of velocity and spatial offsets are shown for three selected 
runs, with green and magenta lines denoting the time spent in off- 
center and disk-crossing quasar phases, respectively. The top pan- 
els show the Model A run with i^k = 440 km s^^, in which the 
BH retains an accretion disk for about half of its wandering time. 
The off-center quasar spends most of its lifetime at low velocities 
(near turnaround) and at spatial offsets of a few hundred parsec. 
The BH crosses the disk only four times, however, so the time 
fraction it spends in disk-crossing quasar phases is very small. The 
runs shown in the middle and bottom panels of Fig. [14] which have 
V < «esc for their respective models, have shorter off-center quasar 
duty cycles and longer wandering times, so the time fractions of the 
off-center quasar phases are very low ~ virtually indiscernible. The 
time fraction spent in disk-crossing quasar phases is more substan- 
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Figure 14. Probability distributions of velocity {left panels) and spatial 
(right panels) offsets for selected runs. In all cases, black lines denote the 
distribution for the entire simulation, green lines denote the distribution of 
time spent in off-center quasar phases only (i.e., t < taccr), and magenta 
lines denote the distribution of time spent in disk crossing phases only. Top 
panels: Model A run with i = 45° and vy^ = 440 km s~^ (trajectory 
shown in Fig.|3j. Middle panels: Model A run with i = 45° and V]^ = 740 
km s~^ (trajectory shown in Fig. W). Bottom panels: Model C run with 
i = 45° and dj^ = 1440 km s~ (trajectory not shown). Note that the 
probability distributions of off-center quasar phases in the higher-velocity 
runs are plotted, but are so low that they are almost invisible. 



tial, but because many disk crossings occur late in the simulation 
as the BH is settling back to the center, these are primarily con- 
strained to the lowest velocity bin (and as expected, the smallest 
distance bin as well). Therefore, while it is certainly possible in 
principle for GW recoil to produce a quasar with a measurable off- 
set, it may prove challenging to observe one directly, especially one 
with a very high kick velocity. 

As mentioned in § [T| |Komossa et al.| ( |2008[ l have proposed 
that SDSSJ0927-I-2943 is a recoiling quasar with a velocity offset 
of 2650 km s~^. We do not plot any recoil trajectories with kick 
velocities this high, because this is well above galactic escape ve- 
locities. Were we to reproduce the distributions in Fig.[T4]for a run 
with Uk — 2650 km s~^, we would find a delta function at Uk, 
because the BH would be ejected from the galaxy with little decel- 
eration. The quasar lifetime, however, would be short. Taking the 
mass and luminosity values estimated by [Komossa et al.| l |2008[ l, 
AfBH ~ 6 X 10** Mq and L ~ O.liEdd, we find that the BH will 
only retain a gas disk of mass ~ 2%Mbh and radius ~ 0.37 pc. 



10 yr. This is much shorter than the 
estimate using 



Komossa et al. 



( 2008| 



Komossa et al.| ( |2008| l is 



which it will accrete within ^ 
upper limit of ~ 10^ yr that 
A^disk,cj ~ Mbh- If the conclusion of 

accurate, then their discovery is extremely serendipitous, because 
it combines two low-probability events. The first is the merger of 
two BHs with spins, mass ratio, and spin alignment configured in 
an unlikely manner to give a kick near the maximum possible GW 
recoil velocity, and the second is the observation of the recoiling 
BH during its relatively short quasar phase. This combined prob- 
ability is difficult to calculate due to its dependence on unknown 



distributions of BH binary parameters, and it is unclear a priori 
whether this rare recoil event is in fact more or less likely to occur 
than a superposition of somewhat unusual quasars. The possible 
role of gaseous outflows is also unclear. Although outflows cannot 
easily explain all the unusual spectral lines, neither can the recoil 
scenario explain all of the narrow lines without outflows or similar 
phenomena being invoked at any rate. Follow-up observations of 
this intriguing source would help to clarify the situation. 

Our finding that high-velocity GW recoil events are unlikely 
to be observed is consistent with the null result of [Bonning et al.| 
(j2007j, who conducted a search for velocity offsets > 800 km 



s ' m quasar spectra. In Fig. 11 we have plotted taccr in all of 
our models, for the range 100 < Uk < 4000 km s"\ The de- 
cline in quasar duty cycle for high velocities can easily be seen, 
and at the highest velocities, all models have taccr ^ lO'^ yr. The 
duty cycle for an off-center quasar with iik = 800 km s~^ is 
6 X 10^ - 5 X 10^ yr for our models in which Vcsc < 800 km 
s~^ (Models A & B), although the BH decelerates over this time. 
The most massive BHs (Models C & D) retain more mass after they 
are ejected and thus have longer duty cycles, but decelerate more 
quickly since they are on bound orbits at this Uk. Duty cycles are 
somewhat longer for these off-center quasars if the accretion lumi- 
nosity is lower, although of course they then will be dimmer. The 
duty cycles for disk-crossing quasars with v ~ 800 km s~^ are 
shorter, however, ~ 10^ yr. In all cases, duty cycles are shorter still 
if one considers higher Uk. Clearly, the probability of observing a 
recoiling quasar is low, although future searches sensitive to smaller 
velocity offsets will increase these probabilities. Finally, as a corol- 
lary to the arguments on the detectability of disk-crossing quasars, 
we note that [King & Dehnen| \2005\ and [Libeskind et al.| ^2006| l 
have argued that some of the brightest ultraluminous X-ray sources 
(ULXs) might be explained as wandering intermediate-mass BHs 
(Mbh ~ 10^ - IO'^Mq) that were ejected from their host galax- 
ies in previous merger events and that begin accreting again upon 
returning to a dense gaseous region. 

There is another possibility for detecting recoil events in ad- 
dition to observing quasar offsets. 'Lippai et al.' f2008}, [Shields &| 
Bonning (2008), and Schnittman & Krolik (20081 have recently 
pointed out that strong flares may occur following a GW recoil 
as the surrounding gas is disrupted and gas marginally bound to 
the BH falls back onto the galactic disk. [Lippai et al.| ( |2008| l ar- 
gue that strong shocks are a natural consequence of perturbations 
to highly supersonic gas, as is the case with thin accretion disks. If 
the shock propagates at Ushock ^ ^k, it will produce a flare with an 

emission spectrum peaking at fc Tshock = (3/16)/i mn I'shock ~ 
19eV (ushock/lOOkms"^)^. Such flares would be observable in 
soft X-rays for recoil kicks > 200 — 300 km s^ ^ ; below these ener- 
gies, most of the emission would be absorbed by interstellar hydro- 
gen. The highest kick velocities ( > 1000 km s~^) could produce 
hard X-ray flares, though these would be short-lived. These photon 
energy ranges would change if the shock propagates at only a small 
fraction of the kick velocity, as Lippai et al. (2008 1 have estimated; 
in this case some fl ares might be observable at UV energies. |Shields| 
& Bonning i 2008 1 estimate a typical flare timescale tflare ~ 10 yr 

1 in 10^ quasars. Schnittman & 



and an observable flare rate of ~ 1 in lU' quasars. 

|Krolik| ( 2008 ) suggest that because the disk is optically thick, recoil 
flares may produce long-lived infrared afterglows, of which up to 
~ 10^ might be visible today. These numbers are uncertain, how- 
ever, and detailed hydrodynamical simulations of recoiling disks 
will be needed to more accurately determine the frequency and du- 
ration of these flares, as well as whether their spectra and lumi- 
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nosities can be distinguished easily from concurrent AGN activity. 
Tlie important point to note, hiowever, is tliat a recoiling BH mak- 
ing numerous passages through the disk may create a succession of 
flares, rather than just an initial flare during the recoil event. This 
could substantially increase the rate at which recoil flares occur. 



3.5 Sensitivity of results to q parameter 

As mentioned above, we have accounted for the ellipticity of the 
gas disk potential $(R) by introducing a small perturbation, g, 
in the radial coordinate such that R = -y/r^ + (z/q)'^ and $ — 
^{R). We use q — 0.99 as a fiducial value for our runs, and here 
we investigate the sensitivity of our results to this choice of q. A pri- 
ori, one can expect the largest perturbation to the BH orbit to occur 
near turnaround, where the velocity is low and the BH is most sen- 
sitive to acceleration from the (perturbed) gravitational potential. 
Therefore, the largest cumulative perturbation should occur for or- 
bits with the most turnarounds, i.e. those with large Uk. Orbits with 
large «k also have longer turnaround times, since turnaround can 
be defined as a fixed fraction of the orbital period, which increases 
for higher kick velocities. Fig. [15] shows the total mass accretion, 
AM/Mbh, as a function of q for low, moderate, and high values of 
V]i. Indeed, we see that very little variation with q occurs except in 
the highest velocity case, < «csc- We can thus conclude that the 
value of q does not matter for simulations with low-to-moderate 
Vk. For the largest kicks, we can assume that choosing q fairly 
close to unity is appropriate, since turnarounds will occur mostly 
at large radii where the degree of ellipticity of equipotential con- 
tours should be small. Furthermore, BH trajectories for Wk ^ v^sc 
are already the most susceptible to other inaccuracies due to phys- 
ical effects that were not included here. For example, at the largest 
radii probed by orbits bound to the baryonic galactic component 
(generally tens of kpc), BHs are likely to encounter perturbations 
due to the triaxiality of the DM halo. The host galaxy is also likely 
to evolve and undergo star formation during the the long ( > Gyr) 
wandering time of the BH in these cases. 



4 DISCUSSION 

We have considered SMBH accretion without requiring that the 
SMBH remains stationary at the center of its host galaxy. Instead, 
the SMBH was assumed to be "kicked" from the center at 100s 
- 1000s of km s~\ which can occur due to GW recoil following 
a BH merger. To calculate the trajectories of recoiling SMBHs in 
a galactic potential, we have developed an analytical model that 
includes a stellar bulge and a multi-component gaseous disk. The 
SMBH mass and gas fraction are free parameters that we varied in 
four different cases (Models A-D). For bound trajectories, we ex- 
amined the dynamics and the wandering timescales of the SMBHs, 
as well as the effect of this wandering on SMBH growth. 

Surprisingly, the kick velocity, i^k, has little effect on the 
amount of mass accreted onto the SMBH while it wanders through 
the galaxy. For all cases, the fractional SMBH growth, AM /Mbu, 
is roughly 10% over the time from the recoil kick to when the BH 
settles back to the galactic center. This is primarily the combined 
result of two different effects: the total wandering time is short for 
low Vk, where the BH quickly settles back to the center, and the ac- 
cretion time is short for high Uk, where the accretion disk bound to 
the BH is small. Due to these competing factors, the timescale for 
off-center accretion, and hence the amount of growth in this time, 
remains almost constant for all Wk- Additionally, bursts of accretion 
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Figure 15. AM/Mbh as a function of the ellipticity parameter q, using 
Model A and i = 45°. The low (-Uk = 100 km s^^, top panel), moderate 
(tik = 440 km s~^, middle panel), and high- velocity (v]^ = 740 km s~^, 
bottom panel) kicks shown in § |3.1| are chosen again here to demonstrate 
variation with q. For clarity, the value log(l — q) is plotted on the x-axis. 
The solid dot denotes our chosen fiducial value q = 0.99. 



during subsequent disk crossings contribute little to the total BH 
mass. The consistency of AM/Mbh implies that GW recoil is an 
effective means of self-regulation for SMBH growth. If the SMBH 
remained stationary at the galactic center, substantial growth could 
occur. This growth could be regulated via more conventional pro- 
cesses such as energy or momentum feedback, but these depend 
on the details of the gas physics and may vary greatly in different 
galactic environments. In contrast, GW recoil depends only on the 
parameters of the merging SMBH progenitors. It is also interest- 
ing to note that because other feedback processes will be halted 
when the recoiling SMBH has depleted its reservoir of bound gas, 
these two methods of growth regulation are complementary. The 
possibility that both wandering and stationary SMBHs self-regulate 
their masses indicates how recoiling BHs could be consistent with 
the observed Mbh — o"* relation, even if large recoils < Wcsc are 
not rare occurrences. The latter conclusion assumes, of course, that 
the escape fraction of recoiling SMBHs is still small; otherwise 
we would expect to see more scatter in the Mbh — cr* relation 
l |Libeskind et al.|2006[ |Volonteri|2007^ or galaxies without a cen- 
tral SMBH in the local universe, which are not observed (e.g., Kor- 



'men dy & Richstone|1995[|Richstone et al.|1998[|Ferrarese & Ford 
2005). 

Even small recoil kicks can have a significant impact on ac- 
cretion and feedback processes. For the lowest kick we consider, 
Wk ~ 100 km s~^, the BH wanders in the galaxy core for afew Myr 
before settling back to the center and carries along a large accre- 
tion disk. Thus, feedback from moving sources will be distributed 
over much larger volumes than the emitting regions of stationary 
AGN. For example, in our fiducial model, a BH kicked at 100 km 
s^^ reaches ~ 30 pc before turning around. Assuming most of the 
AGN emission originates from a region ~ 10^ — 10"^ rs across, 
the initial volume in which energy is released by the BH increases 
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by about 8 orders of magnitude in this case. Most of the energy 
will be deposited at the largest radii, since turnaround takes the 
largest fraction of the orbital time. This implies an even greater dis- 
parity between the impact of moving and stationary BHs on their 
surrounding environments. 

The importance of GW recoil in the evolution of SMBHs 
clearly depends on how often significant recoil kicks occur. This 
is difficult to estimate, because the recoil kick distribution depends 
on the spin and mass ratio distributions of the progenitor BH bina- 
ries, which are not well known. 'Schnittman & Buo nanno| P007| l, 
ICampanelli et al.| l|2007a |l, & Baker et al. (20081 have calculated 
kick distributions; Schnittman &. Buonanno (2007) use an approx- 
imate method, and the latter two derive empirical formulae from 
numerical relativity simulations. The differences in the results of 
the latter two stem from different assumptions of how the kick ve- 
locity scales with BH progenitor mass ratio. Assuming BH spins 
of a = 0.9, mass ratios in the range 1/10 < q < 1, and ran- 
dom spin orientations, the three groups respectively estimate that 
(12, 36, & 23%) of recoiling BHs will have vy, > 500 km s'^ and 
(3, 13, & 9%) > 1000 km s"\ Note that these fractions are not 
insignificant, although they might be much smaller if BH spins are 
on average lower or are preferentially aligned (cf. |Bogdanovic et al.| 
|2007[ ). To constrain the actual recoil kick distribution, better obser- 
vational constraints on merging BHs or observations of individual 
recoiling BHs will be required. 

Our results indicate that the latter method for constraining GW 
recoil, i.e. direct observations of recoiling BHs, will be challeng- 
ing, because several factors make recoiling BHs difficult to observe. 
First, the BH must be actively accreting to appear as an observable 
source. The recoil could be observed either via spatial offsets from 
the galactic nucleus or via velocity offsets in the quasar spectrum. 
In both cases, larger offsets resulting from large recoil kicks are 
more easily observable, but large kicks imply a small disk bound 
to the BH and thus a short off-center quasar duty cycle. After it ex- 
hausts its supply of bound gas, the moving BH could again become 
an active source on subsequent passages through the galactic gas 
disk, an effect which could be observed as a velocity offset in its 
spectrum. This could also be a source of periodic quasar or disk 
flare activity. However, the greater the velocity offset, the shorter 
the duty cycle for these "disk crossing" quasars, since they pass 
rapidly through the relatively thin galactic disk. It is worth empha- 
sizing that because our results indicate that GW recoil may be diffi- 
cult to observe directly and that the growth of wandering BHs may 
be consistent with observed galaxy-BH relations, neither of these 
observational methods may be able to rule out a significant popula- 
tion of recoiling BHs. 

If, indeed, GW recoil is a relatively common phenomenon in 
SMBH mergers, this represents a fundamental change in our un- 
derstanding of SMBH-galaxy coevolution, because we cannot gain 
a complete picture by considering only stationary SMBHs. In re- 
ality, GW recoil and AGN feedback are likely to work together to 
produce the observed SMBH-galaxy relations. 
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